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A Lie-Poisson bracket is presented for a five-field gyrofluid model, thereby showing the model to be Hamilto¬ 
nian. The model includes the effects of magnetic field curvature and describes the evolution of the electron 
and ion gyro-center densities, the parallel component of the ion and electron velocities, and the ion tem¬ 
perature. The quasineutrality property and Ampere’s law determine respectively the electrostatic potential 
and magnetic flux. The Casimir invariants are presented, and shown to be associated to five Lagrangian 
invariants advected by distinct velocity fields. A linear, local study of the model is conducted both with and 
without Landau and diamagnetic resonant damping terms. Stability criteria and dispersion relations for the 
electrostatic and the electromagnetic cases are derived and compared with their analogs for fluid and kinetic 
models. 


I. INTRODUCTION 

Reduced electromagnetic fluid models constitute versa¬ 
tile tools for the study of multi-scale phenomena includ¬ 
ing in particular the interaction of turbulence with mag¬ 
netohydrodynamic perturbations exhibiting meso-scale 
structuresji Examples include magnetic islands^i^ edge 
localized modesj^i^ resonant magnetic perturbations^ 
as well as fishbone^ and Alfven modes.— Irrespective 
of the phenomenon that a particular fluid model aims 
to describe, the underlying system of charged particles 
interacting with electromagnetic fields is a Hamiltonian 
system and in addition to energy and other invariants 
related to symmetry properties, it may possess approx¬ 
imate Poincare or adiabatic invariants such as wave ac¬ 
tions. It is highly desirable that a model giving a reduced 
description should retain this important property in the 
ideal limit. By ideal limit, we mean the limit of the 
model when all dissipative terms, such as collisions, Lan¬ 
dau damping and dissipative anomalous transport terms 
are neglected. Thus, Hamiltonian systems conserve en¬ 
ergy for closed boundary conditions and the Hamiltonian 
formulation is useful for investigating the local properties 
of the dynamics that are independent of the drive. 

Casting a system into its Hamiltonian for m^^d^ confers 
several practical advantages. One of the most important 
is the existence of families of invariants, called Casimir 
invariants, which are found in noncanonical Hamiltonian 
systems due to the degeneracy of the cosymplectic ma¬ 
trix. The functional that results from the addition of 
the Casimirs to the Hamiltonian has non-trivial equi¬ 
librium states as stationary points. In the absence of 
a Poisson bracket, by contrast, the existence of non¬ 
trivial e quil ibrium states is not guaranteed. For exam¬ 
ple, Ref. presents an example of a seemingly reason¬ 
able fluid model that lacks physical equilibria with closed 
streamlines because the equilibrium equations imply that 
some fields are multiple-valued on closed streamlines. 
We can also take advantage of the Hamiltonian formu¬ 
lation to construct “energy principles” for the investiga¬ 
tion of the stability of such non-trivial equilibrium states 
by examining the second variation of the aforementioned 
functional— Another advantage is that imposing con¬ 


straints on a system is straightforward in the Hamilto¬ 
nian formalism— i Lastly, the Hamiltonian formalism can 
be used to facilitate the calculation of the statistical av¬ 
erage of the zonal flow growth rate— ^ 

Among the several classes of fluid models, of particular 
importance are the ones that retain the effects of finite 
ion temperature, principally for describing instabilities 
with growth rates comparable to the ion diamagnetic fre¬ 
quency or modes with perpendicular wavelengths of the 
order of the ion Larmor radius. Whereas “cold ion” mod¬ 
els have been shown to possess noncanonical Hamiltonian 
formulationsthe task of formulating such “hot-ion” 
models that satisfy the Hamiltonian property has proven 
difficult. For example, efforts to identify the Hamilto¬ 
nian structure of the four-field model of Ref. iUw ere un¬ 
successful, even though it conserves energy—^ The main 
difficulty with such models lies in the nonlocality of the 
ion dynamics caused by Larmor gyration. One way to 
approximate nonlocal terms is by a Taylor-series, us¬ 
ing k±pi as a small parameter. An exampfo of such 
a so-called FLR model was given in Ref. 1^, where a 
Hamiltonian four-field model is constructed, using the 
“gyromap” technique to introduce finite ion temperature 
into the cold ion limit of Ref. [U. Unfortunately, we are 
unaware of any numerical implementation of this model, 
possibly because it requires high-order derivatives and, 
consequently, additional boundary conditions. 

An alternative approach for constructing fluid mod¬ 
els with a finite ion temperature is to truncate the mo¬ 
ment hierarchy of the gyrokinetic equation.— This 
leads to the use of nonlocal averaging operators that ac¬ 
count for the full range of perpendicular wavelengths. 
The resulting models are called gyrofluid models. Sur¬ 
prisingly, gyrofluid models are more readily amenable 
to Hamiltonian formulations than FLR models. Exam¬ 
ples of Hamiltonian electromagnetic gyrofluid models are 
given in Ref. [2^ for an incompressible (three fields) and 
Ref. [lO for a compressible (four fields) model. The four- 
field gyrofluid model advances the first two moments of 
the distribution function for each species, or the ion and 
electron densities and parallel momenta. Zacharias et al. 
have shown that simulations of magnetic reconnection 
using this model are in good agreement with gyrokinetic 
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simulationsand Comisso et al. have used it to bring to 
light the acceleration of magnetic reconnection by non¬ 
local gyrofluid effects.— Grasso et ai, by contrast, have 
used it to examine the stabilizing effects of ion diamag¬ 
netic drifts on the growth and saturation of tearing modes 
in inhomogeneous plasmai^ 

In the present paper, we propose a Hamiltonian five- 
field electromagnetic gyrofluid model that is an exten¬ 
sion of the model presented in Ref. M The new model, 
like its predecessor, is a truncation of a more complete 
one proposed by Snyder and Hammett, which advances 
six moments for the ions and two moments for the elec¬ 
tron dynamics— We note that Scott— >2^ has shown that 
achieving energy conservation requires modifying several 
of the terms in Ref. [ 2 ^ involving higher order moments. 
We will likewise show that constructing a Hamiltonian 
model requires modifying the terms involving the higher 
order moments in our model. The new model extends 
that in Ref. [sO by the addition of the evolution of the ion 
temperature. As in the previous model, ion compressibil¬ 
ity effects and field curvature are also included, allowing 
it to describe ITG, KBM, drift waves and tearing modes. 
To demonstrate the properties of the model, we present a 
linear, local study of electrostatic slab ITG and toroidal 
electromagnetic ITG modes. 

The paper is organized as follows: In Section|H]we give 
the normalizations of our variables and present the ideal 
limit of the dynamical model. In Section IHII we give the 
Hamiltonian formulation of the model equations by pro¬ 
viding a conserved energy that serves as the Hamiltonian 
and a Lie-Poisson bracket that satisfies the Jacobi iden¬ 
tity. In Section IIVI we calculate the Casimir invariants 
of our system and from them, in Section |V] we construct 
five “normal fields” which are field variables in which the 
dynamical equations and the bracket take a very sim¬ 
ple form. Lastly, in Section Ivll we perform a local, linear 
study of the model with particular emphasis on the study 
of the ITG and KBM modes. We present stability criteria 
for both the ideal model and a model with linear dissi¬ 
pation terms representing the effects of parallel Landau 
damping and the drift resonance. We investigate several 
well known stabilizing factors of the instability to show 
qualitative agreement with kinetic models. 


II. IDEAL MODEL 

We first present the ideal portion of our model by 
omitting collisional diffusion and wave-particle interac¬ 
tion terms, which will be examined in Sec. ED 

We are interested in a model that describes the desta¬ 
bilization of the drift wave excited by the ion temperature 
gradient. Due to the acoustic nature of the instability, 
we cannot neglect ion motion along the field lines; there¬ 
fore, we keep ion compressibility effects. Also, because 
we want to investigate toroidal plasma with finite /3, we 
include electromagnetic effects. Lastly, to represent the 
influence of toroidicity, we allow for magnetic curvature. 


We consider the evolution of the the magnetic flux tjj, of 
a magnetic field Ti ~ 2 . + Vip x z, the ion density Ui, the 
parallel velocity of the ion guiding centers Ui = z-Vi, the 
electron density rie and parallel velocity Ug = z • Ve, the 
electrostatic potential (j) and the parallel ion temperature 
Ty. We normalize these quantities in the following way: 

{ni,ne,1p,(l),Ui,Ue,T\\) = 

^ ^ j) ecj) ^ ^ ^ 

5 5 5 rji 5 5 7 rf~\ 

'^o Pi^o '^ti '^ti 

where the carets denote the dimensional variables. Here 
rio, Bo and Tg are the background density, magnetic 
field and electron temperature, pi = vulojd is the ion 
Larmor radius, where vu = is the ion ther¬ 

mal speed, ujci = eBojmi is the ion cyclotron frequency, 
Ln = no/|Vn| is the density scale-length and r = TgjTi 
is the ratio of the species temperatures. We also normal¬ 
ize the independent variables according to: 




(t, fc||,/Ci) 


tVti 


, fcll Lrfl , k^P', 


( 2 ) 


With these normalizations, our evolution equations are 
as follows. The equations that describe the ideal evolu¬ 
tion of ion quantities are 


dui d , 

= -V||Ui - 2ud^{ni -I- + T||), 


dt 

d(d' -I- Uj) 
dt 


= —ViiTii — Viini — Aud 


dui 
dy ’ 


^ = -(7 - 1 )V||U* - + $ + T||) 


( 3 ) 

( 4 ) 

( 5 ) 


whereas the equations describing the evolution of electron 
quantities are 


dug 

dt 


-VllMe -b 2-lid—(Ue - (/)) , 
dy 


d{tp - pUg) 

dt 


1 ^ r. 

-V\\ng -b 2p,Ud^:— 
T dy 


( 6 ) 

( 7 ) 


In Eqs. @“(0, df/dt = df/dt+ [<&,/] and Vy/ = 
dfjdz — [4',/], with [•,•] denoting the canonical Pois¬ 
son bracket, so that [f,g] = z • (V/ x Vg). Also, 7 is the 
adiabatic index, Ud = L^/R is the normalized curvature 
drift velocity, R is the radius of curvature of the magnetic 

field and <I> = are the gyro-averaged 

1 /2 

(j) and tj). The symbol Vo refers to the gyroaveraging 
operator introduced in Ref. [U and is defined by 


rn = exp(ivi)/y7-vi){, (8) 


where Iq is a modified Bessel function of the first kind 
and the result of Eq.([5]) should be interpreted in terms of 
its series expansion. At this point, we note that only the 
ion guiding centers respond to the gyroaveraged value of 
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the electromagnetic field. Therefore, we are required to 
use the gyroaveraged value of the electrostatic potential 
in the E x B drift advecting the ions whereas, electrons 
are advected only by the local value of their E x B drift 
since we neglect the electron Larmor radius. 

Equations ©-dll) are closed by the parallel component 
of Ampere’s law 

= = + (9) 

TPe 

with j = z • J being the z-component of the current den¬ 
sity, and by the quasineutrality condition 

He = + (r„ - 1 ) </., ( 10 ) 

with To = . Here, is the gyrophase- 

independent part of the real space ion particle den¬ 
sity and the (Eq — 1 )^ term comes from the gyrophase- 
dependent part of the distribution function. It repre¬ 
sents the ion polarization density due to the variation of 
the electric field around a gyro-orbit. We leave /3e unre¬ 
stricted so that we can describe both “inertial” (/3e “C /.t) 
and “kinetic” {(3e ^ p) Alfven waves. Since our only 
temperature equation involves the parallel temperature, 
from now on we will drop the subscript from . 

It is interesting to compare the model presented in 
equations© -d3 to one obtained from the models of 
Refs. l26l - l^ by discarding all the terms involving high- 
order moments and associated terms. By “associated” 
terms, we mean for example that discarding Tj^ requires 
that one also discard terms involving the gyroaveraging 
operator Ji, since the latter terms result from the ef¬ 
fects on gyroaveraged quantities of the variations in the 
perpendicular temperature. The link between T±^ and 
Ji is reflected in the fact that for energy conservation, Ji 
terms must appear together with Tj^, as noted in Refs.l^ 
andH^ The omission of the terms containing Ji means, 
in effect, that we neglect VJq. Compared to such a trun¬ 
cated model, the Hamiltonian model in Eqs. ©-© lacks 
any tra pped p article effects (terms proportional to VyR 
in Refs. [26l - [2al and has a less accurate treatment of ELR 
terms (due to the omission of the Ji terms). The two 
models also differ in the coefficients of the various curva¬ 
ture terms. In the continuity equation, for example, the 
argument of the curvature operator in the truncated ver¬ 
sion of the model of Refs. [26l - l^ is $ -I- P||/2, while that 
in our model is <& -|- py. This difference is necessary in 
order for the five-field model to conserve energy. In fact, 
we note that the curvature terms in Eqs. ©, © and 
are the same as the ones found in the correspond¬ 
ing equations of the FLR fluid model of Ref. [s^, which 
evolves three ion moments, as we do, and conserves en¬ 
ergy. Lastly, we note that the factor of four in front of the 
curvature term in the momentum equation, Eq. ®, does 
match the corresponding term in Refs. despite the 

fact that for the four-field model of Ref. [^, satisfying the 
Jacobi identity required halving this factor. The conclu¬ 
sion of these observations is that constructing Hamilto¬ 
nian models requires modifying the truncated moment 


expansions, but that the correct terms are recovered as 
on increases the order of the model. 


III. THE HAMILTONIAN FORM 

The system described in Sec. HIl conserves the following 
energy: 

H = ^ (— +nl + + uf 

2 Jp \T 7-1 

+7|;|VV’P + , (11) 

where V denotes the spatial domain of interest and the 
boundary conditions are such that surface terms vanish. 
The successive terms of the functional of Eq. (HU repre¬ 
sent, respectively, the electron and (two terms) ion ther¬ 
mal energies, the parallel component of the electron and 
ion kinetic energies, the magnetic energy and the electro¬ 
static energies of ions and electrons. Taking the energy 
functional as the Hamiltonian of our 5-field model, we can 
write the set of equations in a noncanonicaU^ Hamilto¬ 
nian form 

^ = {e,iL}, * = !,...,5, (12) 

with being the field variables and {•,•} being a non- 
canonical Poisson bracket. We employ the dynamical 
variables rii, Mi,ne, Me,T, where Mi = + Ui is 

the canonical ion momentum and Mg = ip — pue, the 
electron one. Additionally, we define fii = rn — 2udX, 
he = Ue — 2udX and T = T — 2udX for convenience. In 
these variables, the bracket given by 

G} = J d^x(^ — fii(, GjjJ -|- [FMi,GMi] 

+ [Ft , Gt \) — Mi ( [Fmi , Gm ] + [An., ] 

-I- ([At, GmJ + [Am^, Gt])) 

— A([A„^, Gt] + [At, G„J -I- [Fmi, Gm*]) 

+ fie([A„^, G„J -I- p[Fm^,GmS) 

+ Me{[FM^,Gn^] + [A„^, GmsI) 

— {FMidzGrii — GMidzFm) 

— {FxdzGMi — GrdzFMi) 

+ {FMedzGrie — GMedzFnS'j ( 13 ) 

satisfies the formulation of Eq. m for the Eqs. ©-©, 
is bilinear, antisymmetric and satisfies the Jacobi iden¬ 
tity. In the above bracket, we have taken 7 = 2 because 
this is the only value of the adiabatic index that allows 
the bracket to satisfy the Jacobi identity, as shown by a 
direct pro of of the Jacobi identity using the techinques 
of Ref. [Il|. The Jacobi for this case will become evident 
in Sec. |Vl 
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IV. CASIMIR INVARIANTS 

One of the most important properties of noncanonical 
Hamiltonian systems is the existence of Casimir invari¬ 
ants, that is, constants of motion for any choice of Hamil¬ 
tonian. A Casimir invariant C thus needs to satisfy the 
relation {F,C} = 0 for any field F. Here, we will set 
dz = 0. The generalization is straightforward. 

Assuming a Casimir functional C{ni, Mi,T,ne, Me) 
and applying the condition {^j,C} = 0 with = Ui, 
^2 = Mi, ^3 = T, ^4 = rie, Cs = Me gives following: 

[n, - 2udX, + [M„ CmJ + [T - 2udX, Ct] = 0 (14) 
[ui - 2udX, CmJ + [M^, Cm] 

+ [T- 2udx, Cm,] + [M„ Ct] =0 (15) 

[m - 2udx, Ct] + ]T- 2udx, Cn,] + [Mi, Cm,] = 0 (16) 
K - 2udx, Cm] + [Me, Cm,] = 0 (17) 

^i[ne-2udX,CM,] + [Me,Cm]=Q ■ (18) 

For the rest of this section, we employ the previously 
defined variables hi, he, T . In addition, we observe that 
F^ = F^. From (I14|) and (IT51) we retrieve no information 
since they are automatically satisfied for any choice of C. 
However, from (fT5)l we get 


g) = 0 and by application of the method of charac¬ 
teristics, we recover the other characteristic direction, 
C = (g{f + hi,Mi) + f{hi-f)Y Finally, employing 

(1^ we arrive at the wave equation d\j.g — 2d‘^mT9 — 0- 
Invoking the method of characteristics once more, we de¬ 
rive the following general form for the Casimir invariants 
corresponding to the ion piece of the bracket: 

C,= J g±{f + h,± V2M,) + f{h, - f). ( 27 ) 

For the Casimir invariants that correspond to the elec¬ 
tron part of the bracket, we need only solve (E5)) to obtain 

Ce = j h±{Me ± ,/JIhe). ( 28 ) 

Thus, a general family of Casimir invariants is given 

by 

C{n,,M,,T,ne,Me) = J g±{f + hi ± \/2Mi) 

+ f{hi-f) + h±{Me± y/Jlhe), ( 29 ) 
where g±, f and h± are arbitrary functions. 


[hi, Mi][CMiM, Cmn, Ctti,) 


+ [Mi,f]{CniT — CMiM, + Ctt) 

+ [h„f]{CM^T-CM,,n,) = Q, (19) 

from (HU) we get 

[hi,f]{CTT — Crum) 

F[f,Mi]{Cm,Mi — Cmit) 

+ [hi, Mi]{CTMi — CMifii) = 0, (20) 

and from (HU we get 

[he, Me]{flCM,M, " C^n,) = 0 . (21) 

Accordingly, we have the following set of equations: 

CMiM, — Cmm — Ctu, = 0 ( 22 ) 

CMiM, — Cn,T — Ctt = 0 (23) 

Ctm, — CMiu, = 0 (24) 

Ctt - Cmn, = 0 (25) 

^CmsM, Cn,rie ~ 9 , (26) 


V. NORMAL FIELDS 

The general form of the Casimir (|29)l suggests the in¬ 
troduction of a new set of variables which are called “nor¬ 
mal fields” (see e.g. Refs. IsM andls^: 


Vi,± = f + hi± V2Mr 

(30) 

Vij =hi - f 

(31) 

Ve,± =Me ± ^JJihe . 

(32) 


We claim that if we express the equations of motion ([3]) 
- © and the bracket of (US in terms of these fields, 
they will take a simple form. To do so, the following 
chain rule expressons for functional derivatives in terms 
of these new fields are required: 


Fm 

=Fv,.+ + 

Fv,, 

/ + Fv,,- 

(33) 

GTii 

=Fvi,+ + 

Fv., 

1 

1 

(34) 

Fm, 

=V2 {Fv, 

,+ ~ 

Fv,,-) 

(35) 

Fm, 

=Fv,,+ + 

Fv, 

_ 

(36) 

Fn, 

=\fh {Fv, 


-Fv„-) . 

(37) 


which must be satisfied by any Casimir invariant. 

We start from Ea. (0^ and integrate it w.r.t Mi to find 
Cm = Ct + f{hi,T). By using the method of char¬ 
acteristics on this result, we infer that the solution has 
the form C = (^g{T + hi,PIi) + f{hi,f)'^, where the () 

symbol implies an integral over the volume of interest. 
Subsequently, we substitute this form of the Casimir into 
(051) to obtain the wave equation d^.{f + g) — + 


Using (I551) - (I57)) the Poisson bracket of (IT^ becomes 

{F,G}=- 2 (V.,/[Fv.,,,Gv,,,] 

+ 2 (V,,+ [Gv.,+ ,Gv,,J + V,,_[Fv.,_,Gv.,_]) 

- {Ve,+ [Fv„^,Cv,J - Ve,-[Fv„..Gv,J) 

+ 2V2 {Fv,_^dzGv,,+ - Fv,__dzGv,,_) 

- y/Jl (Gv„,+ i9zGv^,+ - Fv,__dzGv,,_)) ■ (38) 
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This simple form of the bracket is called a direct 
product^, and its form immediately ensures the Jacobi 
identity. Since the inner brackets satisfy the Jacobi iden¬ 
tity, so do their sums which constitute the larger bracket 
of Eq.dni). 

Having expressed the bracket in terms of the normal 
fields, we can now write down the equations of motion 
that these fields satisfy, viz. 

+ [A.±, v,.±] ± V2d,A,± = 0 ( 39 ) 

+ [Ae,± , Ve,±] T =0 (40) 

^ + [A,/,V,y]=0, (41) 

where 

Ai,± =$ -I- rii -I- T ± (42) 

A,j =<^ + ni-T (43) 

Ae,± = ± Ue (44) 


Parameters x ^-nd i' of the added dissipative terms are 
tuned so that the response function of a gyrofluid model 
matches the kinetic one in the slab and the toroidal lim¬ 
its, respectively. Their values have been computed in 
Refs. ^ andl^ and found to be y = ^ and v = 2.019. 
Although the x value is exact, the numerical value of 
v has not been calculated for the particular model we 
are presenting but for a similar gyrofluid model. Nev¬ 
ertheless, we will adopt it. The reason is that here, we 
are mainly concerned with the non-dissipative, Hamilto¬ 
nian part of the model and the addition of the dissipative 
terms is not intended to enhance the accuracy of the re¬ 
sults, but merely to show the reader that such a modifica¬ 
tion is indeed possible. Correct treatment of dissipation 
would require the proper study of the response function 
of a kinetic model containing the same physics and the 
numerical minimization of the error in matching it with 
the response function obtained by ©-([T]). Such a study 
is beyond the goals of this paper. 

The linearization of the equations of motion and the 
closure relations in Fourier space result in the following 
system of equations: 


are stream-functions that simply convect the fields 
Vs,±//. The latter are therefore Lagrangian conserved 
quantities. Note that in a turbulent system, equiparti- 
tion results in the flattening of the profiles of Lagrangian 
invariants!^ 


1 /2 '' 

-Lofii =w*rQ' (&)</) - kyBoyUi 

— -|- Tq^ (^)'^ T 

- kzUi , (45) 


VI. LINEAR STUDY 

In this section we linearize ©-O and the two closure 
relations ©-(nni) about an inhomogeneous equilibrium 
configuration. Then, after deriving the dispersion rela¬ 
tion, we study the linear stability of the ITG mode. We 
assume that the densities and temperature vary linearly 
in the x direction, i.e., that these quantities have the form 
/ = x/Lf + Sf with Sf = /exp(ik-x —iwt). This may be 
interpreted as a local study, in the WKB sense, for modes 
satisfying k±L± ^ 1 and A:||L|| ^ 1, where L± and Ly 
represent equilibrium scale-lengths. Our purpose is to 
obtain some physical understanding of our model and 
see how accurately it can describe the various modes of 
interest. Next, we assume ^eq = 0 and Vipeq x z = Boyij 
with Boy = 8 - constant and Ui^eq = 0. We note that 

in Fourier space, the operator To is ro( 6 ) = e~^Io{b), 
where b = (or 5 = in our normalized units). 

Fven though we mentioned that the model is Hamilto¬ 
nian only for the choice 7 = 2 , in the following we keep 
7 general to investigate its effect on the behavior of the 
modes and we subsequently set 7 = 2, to recover the 
results for our model. 

Moreover, we add two dissipative terms to Eq. © 
that are related to the parallel and toroidal resonances. 
Therefore, from now on, we make the distinction between 
the non-dissipative, i.e. Hamiltonian, gyrofluid model 
and the one where dissipation terms are included. 


-u;(Tl^^{b)'tp + Ui) = - kyBoyf - uj^r]iTl^^{b)tjj 

kyBoy^Q (^)0 kyBoyTli 

— k^fii — k^T 

-k^ryyb)^, (46) 

-ujf =u}^r]iTy'^ {b)^ - ( 7 - l)kyBoyUi 

-2uj,€{n^+ryyb)^+f) 

— {-j — l)kzUi+ 2w\uj^\ef (47) 
+ ix\k\\\f, 

UJTIq — kyBoyUo 

+ 2a;,e ^ - k^iie , (48) 

~ - ky 

-Uj{lp - pUe) = - kyBoy4> H-V' J- Boyfle 

T T 

^ hz ^ 

-I- 2a;*epUe — k^cj) H- fie , (49) 

T 

fie =ry^( 6 )hi -I- (ro( 6 ) - 1)$ , (50) 

= -Ue + f'y^{b)ui. (51) 

Pe 

Note that {h)Boy = Boy and, to be clear, recall the 
ion and electron density and parallel temperature gradi¬ 
ents vary linearly, i.e., = x/Ln^^Ue = and T = 
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x/Lt- We simplify the result by setting + Boyky 

and by defining the parameters r]i = LnjLT-, e = UdL^, 
and TYi — Also, — i^cT^jcBo^ikyJLji^ is the 

usual diamagnetic frequency. In dimensionless variables 
it is expressed as w* = ruuky/Ln- 

A. ELECTROSTATIC DISPERSION RELATION 

The electrostatic limit, which is applicable for low-/3 
conditions^ leads to a cubic dispersion relation that of¬ 
fers the opportunity of comparing analytic solutions of 
the gyrofluid model to kinetic results. To make contact 
with well-known analytic results for the slab branch of 
the ITG mode, we also neglect toroidal effects. That is, 
we drop all toroidal terms of Eqs. SU-dlgl), set ifj = 0, 
and study the slab, electrostatic ITG modes, where the 
drive is due to the coupling of the parallel transit of par¬ 
ticles with the temperature gradient. We notice that in 
this case, the electron and ion fields are decoupled so we 
only use the ion field of Eqs. (H5]) - (H7)) . along with the 
quasineutrality condition of (ISTl) and the electron adia¬ 
batic response Ug ~ 4^/t. After straightforward manip¬ 
ulations, we obtain a dispersion relation with real part 
given by 

-f 1 - Toib)^ + y/cy - 

- ro(&)w,a;^ -I- ro( 6 )fcj|U;, ((7 - 1) - 77 ^) = 0 (52) 

and imaginary part by 


Returning to Eq. (|5^ we can infer two stability crite¬ 
ria. The first one comes from neglecting the dissipative 
terms, hence having just the real part of the dispersion 
relation and by demanding the third-order polynomial to 
have only real roots. This is done by setting the cubic dis¬ 
criminant equal to zero and by that deriving a quadratic 
equation in rji. To investigate the case of finite k±, we 
obtain the stability criterion by making no approxima¬ 
tion on ro(6). The result is shown in Eig. [T] where rjcrit 
(the root of the quadratic equation mentioned above) has 
been plotted as a function of b for various values of fey. 
The curves depicted in Fig. [1] are qualitatively similar to 


^crit 



FIG. I: Stability criterion with finite k± as given by b. 
Here fcy ranges from 0.05 to I.O 


l + (l-r„(6))r 2 


to “h YQ(b')uj^uj 

kf^il + r)' 


= 0. (53) 


A simple picture of the dynamics of the ITG insta¬ 
bility, as determined by e.g. (15^ . is given in Ref. |4l|, a 
picture that will be helpful for interpreting our results. 
A basic scenario for the development of the instability 
can start with a density perturbation, which is confined 
to variation along the field line because the E x B ve¬ 
locity across the field lines is incompressible. Electrons 
respond adiabatically to this ion density perturbation in 
order to maintain quasineutrality and in doing so set up 
an electrostatic potential. This potential perturbation 
then leads to an A x i? drift that injects cool ions into the 
compressed (increased density) region, thereby lowering 
the pressure. That is, the plasma exhibits negative com¬ 
pressibility. The resulting lowered pressure then draws 
ions along the field lines by generating a uy. Since the 
whole picture develops in time and moves perpendicular 
to both the magnetic field and the temperature gradient, 
the ions that move parallel to the field line, prompted by 
the lowered pressure, end up increasing the initial density 
perturbation. 


those reported in Ref. where a kinetic model was used. 

The second stability criterion we deduce, concerns the 
case of perturbations with very long parallel wavelengths 
and comes from setting the imaginary part of the dis¬ 
persion relation equal to zero, solving for uj under the 
condition fcy = 0, and eliminating it from Eq. (l52l) . With 
this procedure we find 

VcrU = 7 - 1 • (54) 


Observe, the critical value depends on the adiabatic in¬ 
dex. The kinetic result for this limiting case is provided 
in Ref. and is given by 


Vcrit^ = 


1 + 25(1-^ 




(55) 


with Ii{b) and Io{b) being modified Bessel functions of 
the first kind. Note that the adiabatic index in the exact 
moment equation for the evolution of the parallel tem¬ 
perature is 3. In Fig. [2] we plot this relation and the 
corresponding fluid approximation of it and we notice 
that our gyrofluid model has the correct asymptotic be¬ 
havior for perturbations with very small perpendicular 
wavelengths provided 7 = 2. However, had we chosen 
7 = 3, we would have gotten the correct asymptotic be¬ 
havior for very large perpendicular wavelengths, at the 
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cost of a non-Hamiltonian model. Moreover, the choice 
7 = 5/3 gives ry^rit = 2/3, the result for the fluid model 
of Ref. m 

The reason behind this discrepancy stems from the fact 
that our model lacks an equation for the evolution of the 
perpendicular temperature. Therefore, all assumptions 
about the correlation of Tj_ and Tj| are made by the choice 
of 7 (with 7 = 3 meaning T± and Ty are uncorrelated and 
7 = 5/3 meaning Tl = Ty) and remain fixed throughout 
the dynamics. Despite this obvious inflexibility of the 
gyrofluid model, it is evident from Fig. [2] that it still 
remains far superior compared to its FLR counterpart. 



FIG. 2: Comparison of critical rj between kinetic, fluid 
and gyrofluid results for the case A:y = 0 

It is helpful to study the ‘fluid’ limit of Eq. (l5^ . which 
is obtained by setting ro(5) = 1 corresponding to 6 = 0. 
This is the limit of very long perpendicular wavelengths 
compared to the gyroradius. Figure [3] shows the stability 
criterion in this fluid limit for three different values of 
7 , results that were previously obtained in Ref. for 
7 = 5/3, where a heuristic explanation was for given for 
the rjcrit limiting value for very long fcy. 


^crit 



FIG. 3: Stability criterion at the ‘fluid’ limit with r = 1 
for different values of the adiabatic index 

To conclude with the electrostatic slab case, we inves¬ 
tigated the growth rate as a function of r. The condition 


T < 1 or, in other words, Ti > is a well-known stabiliz¬ 
ing factor for ITG, which is of particular importance for 
the hot-ion cores of tokamaks.— Indeed, the behavior 
we found was the expected one. 


B. ELECTROMAGNETIC DISPERSION RELATION 

To be applicable to the higher plasma pressure 
achieved by auxiliary or alpha-particle heating, the the¬ 
ory must include the electromagnetic effect. In fact, this 
effect becomes important at surprisingly low-/? because of 
other small parameters in the problem. It is well known 
that increasing (3 stabilizes ITG modest, but leads to 
the onset of kinetic ballooning modes, also known as the 
Alfvenic ITG modes (AITG)4^ For toroidal ITG modes, 
the drive comes from the coupling of curvature and VR- 
drift terms with the temperature gradient, so that we 
must now keep the toroidal curvature terms. It can be 
easily seen that, to lowest order, the electromagnetic ef¬ 
fect is stabilizing. The electromagnetic perturbation cre¬ 
ates a small component of B that is perpendicular to 
both the background magnetic field and the pressure gra¬ 
dient. This component then leads to the development of 
a force on the ions, parallel to the field lines that opposes 
the attraction from the pressure lowering of ITG. 

In the remainder of this section, we follow the analysis 
of Kim, Horton and Dong^ and compare our gyrofluid 
results with their local kinetic ones. We note, however, 
that complete agreement cannot be expected since Kim et 
al. has one extra parameter, namely rjf.. We also note that 
the eigenfrequencies for the model in Ref. [2^ lie within 
a few percent of the kinetic results, so that comparing 
our model to the kinetic results is effectively equivalent 
to comparing it to the Snyder and Hammett model. 

Because the dispersion relation becomes unwieldy and 
doesn’t provide much physical insight, we refrain from 
displaying it here. Instead, we solve it numerically 
and present the results. Figure 2] shows the normal¬ 
ized growth rates for (a) the ideal and (b) the “Lan¬ 
dau” versions of our model as a function of (3 when 
rii = 2.5, b = 0.5,fc|| = 0.1, e = 0.2, = 1 and t = 1. 

We also provide the kinetic and fluid model results from 
Ref. [43 for comparison. By the “Landau” version we 
mean of course the Hamiltonian model augmented by 
dissipative terms modeling the damping caused by the 
wave-particle interactions. From Fig. Hal it becomes im¬ 
mediately clear that the nonlocal treatment of the ion re¬ 
sponse in the gyrofluid model reproduces the main qual¬ 
itative features of the kinetic result much better than 
the fluid model. Compared to the fluid result, the gy¬ 
rofluid one gives stronger stabilization of the ITG modes 
and lower thresholds for the excitation of KBMs. This 
is related to the toroidal resonance. In both Fig. Hal and 
Fig. |4b]we observe the close connection between the sta¬ 
bilization of the ITG mode and the excitation of the Ki¬ 
netic Ballooning mode in accordance with what kinetic 
theory predicts. The addition of dissipative terms makes 












the curves shift closer to the kinetic result although we 
remark that at low growth rates the agreement is less 
satisfactory. 



(a) Hamiltonian model without dissipative terms and 
comparison with kinetic and fluid results. 

r 



(b) Model with dissipative terms and comparison with kinetic 
and fluid results. 

FIG. 4: Normalized growth rate vs. /3 for 
rii = 2.5, b = 0.5, fey = 0.1, e = 0.2, rn = 1, and r = 1. 
Toroidal ITG is the black line while KBM is the red. 

Here, we pause to explain an interesting effect, the 
destabilization due to the addition of dissipation of two 
previously marginally stable modes (7 = 0). For exam¬ 
ple, for the GF model it is seen in Fig. 0^ that without 
dissipation the KBM becomes unstable at /3 « 0.010, 
while in Fig. 00 it is seen for the same case with dissipa¬ 
tion that this mode is destabilized for all values of /3. A 
similar shift from stability to instability can be observed 
upon comparing these figures for the ITG mode, which 
is seen in Fig. 0^ to transition to instability at a some¬ 
what smaller value of /3. To understand these transitions 
consider Fig. 01 where we plot the real parts of the fre¬ 
quencies versus /3 for four modes of the GF model without 
dissipation. In this figure the two upper most modes are 



FIG. 5: Real frequency vs. /3. All parameters are the 
same as in Fig.4 


marginally stable for small values of (3, then as P ap¬ 
proaches the transition value near 0.010 they collide and 
produce instability. This bifurcation, which is standard 
in Hamiltonian systems, is called the Hamiltonian Hopf 
(or Krein) bifurcatio n^^ds . Observe, the same bifurca¬ 
tion occurs when two marginally stable modes collide as 
P is decreased to a value near to the KBM transition but 
closer to 0.009, producing the unstable ITG mode. (Af¬ 
ter the transitions there are also damped modes that are 
not shown in the figures.) 

The dissipative destabilization observed in Fig. 0b] is 
a generic feature of Hamiltonians systems with negative 
energy modes (NEMs). Indeed, in Hamiltonian systems, 
whenever a Hamiltonian Hopf bifurcation occurs, one of 
the modes must be a NEM, and such modes have the 
property of getting destabilized with the addition of dis¬ 
sipation (see e.g. Ref. [13 for a Hamiltonian version of 
the classical Kelvin-Tait theorent^) . One could perform 
a calculation like those of Refs. Ho and [s^, where the 
modal eigenvector is inserted into the perturbation en¬ 
ergy in order to show explicitly that it is an NEM and 
then show the dissipative terms remove energy from this 
mode, but such a calculation is outside the scope of the 
present paper. (A similar situation happens when energy 
is added to a positive energy mode.) Also note, the pre¬ 
viously unstable ITG and KBM modes of Fig. 0a1become 
less unstable at the onset of dissipation, as is shown in 
Fig.0b|due to the fact that it becomes harder for a mode 
to grow when there is less available energy in the system, 
which is consistent with this scenario. 

We reiterate that the purpose of our model is to im¬ 
prove the nonlinear fidelity of fluid models. From that 
perspective, we view the quality of agreement in Fig. 0| 
as adequate. 

In Fig. 01 we display the dependence of the growth rates 
of the ITG and KBM modes on fey for various values of 
P for the model augmented with the dissipative terms, 
with all other parameters remaining the same as in Fig.01 
Values of P are in the range 0.000 — 0.012. We observe 
that the stabilization through the electromagnetic effect 
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becomes more efficient with decreasing fc||. Further, we 
see again the near simultaneous stabilization of ITG and 
destabilization of KBM as was noted above. 



FIG. 6 : Growth rates of the ITG-KBM modes as a 
function of the parallel wavenumber for the gyrofluid 
model with dissipative terms. 


For large values of fey the mode is stabilized by the 
large parallel ion transit term4^i^i^ Intuitively, we can 
understand that the mode is limited by the fact that an 
appreciable initial density perturbation cannot be cre¬ 
ated within an arbitrarily small length scale. Even before 
this limit is reached, though, the negative compressibil¬ 
ity mentioned above is proportional to the ratio of the 
ion diamagnetic to the sound frequency (wpi/k^^Cs), so 
that coupling to the sound wave acts as a source of sta¬ 
bilization. On the other hand, the initial density and 
potential perturbations, as well as the resulting pressure 
lowering, are all proportional to /cy. Therefore, a finite 
fcy is needed to overcome the stabilizing effect of curva¬ 
ture and /3. Thus, the mode becomes most unstable at 
some intermediate value. We remark that in practice a 
complete treatment of the effect of fcy requires a nonlocal 
approach since fcyTy ~ 1 normally applies, so that the 
WKB approach in insufficient. 

In Fig. [7] we illustrate the behavior of the growth rate 
versus fc_L for various values of /3, with the same parame¬ 
ters as those of the previous figures. We notice that the 
peak growth rate occurs around fc_L ~ 0.65 and does not 
change much with /3. Furthermore, the stabilizing effect 
of /3 is almost uniform for wavenumber values higher than 
this. This could be attributed to a very high phase veloc¬ 
ity of the wave, which leaves few particles with the right 
thermal speed to resonate with it. For smaller wavenum¬ 
bers, however, the stabilization due to /3 becomes ineffec¬ 
tive. This is because in this region the parallel ion tran¬ 
sit term becomes significantly larger than the curvature 
term and becomes the dominant stabilizing effect. An¬ 
other important stabilizing effect at high k^_ comes from 
FLR physics. Namely, the ions respond to the gyroav- 
eraged electrostatic field, thereby reducing the effective 
E X B velocity. 

In Figs. [5] and IHl we compare the Hamiltonian model 
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FIG. 7: Growth rates of the ITG-KBM modes as a 
function of the perpendicular wavenumber for the 
gyrofluid model with dissipative terms. 


augmented by dissipative terms and the kinetic result 
from Ref. Both figures suggest some common fea¬ 
tures: again, the qualitative similarity between the 
Hamiltonian and the kinetic curves is evident. However, 
there is a quantitative disparity since the Hamiltonian re¬ 
sult is roughly three times higher than the kinetic one at 
the peak value of 7 . This deviation seems to be corrected 
by taking into account the dissipative terms which low¬ 
ers the results to at most 30% off from the kinetic ones 
at peak growth rate. This amendment, though, doesn’t 
come without its own problems, namely the erratic be¬ 
havior of the dissipative model at low values of 7 . 



FIG. 8 : Comparison between growth rates of the ITG 
mode as a function of the parallel wavenumber for the 
ideal Hamiltonian model, the model with linear 
wave-particle (Landau) damping, and the kinetic 
results. 


VII. CONCLUSIONS 

We presented a Hamiltonian, five-field, electromag¬ 
netic gyrofluid model that evolves three moments for 
the ions (density, parallel momentum, and parallel tern- 
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FIG. 9: Comparison between growth rates of the ITG 
mode as a function of the perpendicular wavenumber 
for the ideal Hamiltonian model, the model with linear 
wave-particle (Landau) damping, and the kinetic 
results. 


perature) and two moments for the electrons. We gave 
the Hamiltonian formulation of the model by providing 
a suitable Hamiltonian and a Lie-Poisson bracket that 
satisfies the Jacobi identity. For this system, we found 
the families of Casimir invariants and, from them, we de¬ 
fined five normal fields in terms of which the equations 
of motion and the bracket take their simplest form. 

To evaluate the physical fidelity of the model, we per¬ 
formed a local, linear study of the dispersion relation. 
We began by comparing the electrostatic dispersion re¬ 
lation with known analytic results for a slab ITG mode. 
We found that the critical iji for the onset of instabil¬ 
ity is unity for all k±. This is very close to the kinetic 
result for perpendicular wavelengths comparable to and 
greater than the ion gyroradius, but it is only half of the 
exact value (two) for long wavelengths. Ordinary fluid 
models, by contrast, yield good agreement at long wave¬ 
length but predict negative values of 77i,crit at moderate 
and short wavelengths. By ordinary fluid models, we re¬ 
fer here to those that are derived from the Braginskii 
model and other long wavelength expansion procedures. 

We subsequently examined the electromagnetic prop¬ 
erties of the model, including toroidal curvature, by com¬ 
paring the dependence of the growth rate on well known 
stabilizing factors and comparing them with the local ki¬ 
netic result of Ref. We found good qualitative agree¬ 
ment although the two models cannot be directly com¬ 
pared since the latter includes the effects of the electron 
temperature gradient. We leave for future work the task 
of including in the model an electron temperature evolu¬ 
tion equation in a manner that preserves the Hamiltonian 
character. 

Given the wide availability of several high-quality gy- 
rokinetic (GK) codes that have been verified and vali¬ 
dated in a broad array of contexts, it is appropriate to 
reflect on the value of gyrofluid (GF) models. Due to 
its nature as a truncated moment expansions of the GK 
model, a GF model such as the one presented here cannot 
aspire to compete with the latter in any but three do¬ 


mains: speed, ease of use, and by virtue of the first two, 
ability to generate physical insight. The success of the 
TGLF code^^— demonstrates that there is a strong de¬ 
mand for an agile quasilinear GF code to understand and 
interpret experimental observations of turbulent trans¬ 
port. The motivation for the development of the Hamil¬ 
tonian GF model presented here is similar but different: 
it is to provide an equally agile tool to investigate multi¬ 
scale nonlinear problems such as those listed in the in¬ 
troduction. In this context, the linear accuracy of the 
model is of secondary importance compared to assuring 
the proper conservation laws and providing a qualita¬ 
tively correct picture of the nonlinear energy transfers. 
It is worth noting, in this context, that a Poisson bracket 
for a gyrokinetic model, demonstrating its Hamiltonian 
nature, has only recently been constructed^ using the 
newly developed technique of gauge-free lifting.— 

Uncertainty quantification (UQ) provides another ap¬ 
plication for reduced models that is worth mentioning. 
The large number of inputs to gyrokinetic codes make 
comprehensive UQ impractical, but the existence of a re¬ 
duced model opens up new avenues for charting model 
sensitivities and subsequently using Bayesian inference 
with a smaller number of runs of the GK code to selec¬ 
tively refine the predictions and reduce the error bars. 
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